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Abstract 

In its customary formulation for one-component fluids, the Hierar- 
chical Reference Theory yields a quasilinear partial differential equation 
(pde) for an auxiliary quantity / that can be solved even arbitrarily close 
to the critical point, reproduces non-trivial scaling laws at the critical sin- 
gularity, and directly locates the binodal without the need for a Maxwell 
construction. In the present contribution we present a systematic explo- 
ration of the possible types of behavior of the pde for thermodynamic 
states of diverging isothermal compressibility kt as the renormalization 
group theoretical momentum cutoff approaches zero. By purely analyt- 
ical means we identify three classes of asymptotic solutions compatible 
with infinite kt, characterized by uniform or slowly varying bounds on 
the curvature of /, by monotonicity of the build-up of diverging kt, and 
by stiffness of the pde in part of its domain, respectively. These scenarios 
are analzyed and discussed with respect to their numerical properties. A 
seeming contradiction between two of these alternatives and an asymp- 
totic solution derived earlier [Parola et al., Phys. Rev. E 48, 3321 (1993)] 
is easily resolved. 

Keywords: liquid- vapor transitions, non-linear partial differential equations, nu- 
merical analysis, finite differences, stiffness. 
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1 Introduction 



Reconciling the vastly different approaches to fluid structure and thermodynam- 
ics afforded by classical integral equation (ie) formalisms and renormalization 
group (rg) theory, the Hierarchical Reference Theory (hrt, [1, 2, 3, 4, 5, 6, 7]) 
presents itself as a particularly effective instrument for studying the critical re- 
gion and liquid-gas phase equilibrium in simple one-component fluids: For sub- 
critical temperatures, T < Tc, the usual formulation of the theory [1, 7] yields 
density intervals of rigorously flat free energy and infinite isothermal compress- 
ibility kt the boundaries of which are readily identified with the densities and 
Qi of the coexisting gas and liquid phases. The binodal so found terminates at 
some temperature T = Tc and density g = Qc m & liquid-gas critical point char- 
acterized by non-classical, partly Ising-like exponents [2]. And far away from 
the coexistence region of the phase diagram, hrt reduces to one of the standard 
approximations of liquid state theory, viz., the popular scheme commonly known 
under the names of Lowest-Order 7 Ordered Approximation (loga, [8, 9]) and 
Optimized Random-Phase Approximation (orpa, [10]). hrt's unified treatment 
of thermodynamic states so diverse should be contrasted with the limitations 
inherent in approaches based on lEs alone: Close to the critical point these 
generally do not have a solution or else develop major deficiencies, and the bin- 
odal is accessible only by way of a Maxwell construction [11]. RG calculations 
alone, on the other hand, while invaluable for illuminating the scaling relations 
valid asymptotically close to the critical point, generally do not allow one to 
determine non-universal quantities such as, e. g., the loci of the critical point 
and the binodal. A theory like HRT that provides comprehensive structural and 
thermodynamic information both close to and away from the critical point and 
exhibits a non-trivial scaling limit is clearly attractive for applications calling 
for high-resolution data on the behavior of a fluid in the critical region. 

In the present two-part series of reports we want to have a closer look at the 
solution of the hrt equations for thermodynamic states of diverging compress- 
ibility, i. e., at the critical point and at phase coexistence. Our motivation for 
this inquiry is twofold: First of all, we aim to extend our understanding of the 
way in which hrt achieves its remarkable description of criticality and phase 
separation beyond mere invocation of its conceptual ingredients, RG theory and 
thermodynamic consistency {v. i.) in particular. Instead, it is on the level of the 
partial differential equation (pde) itself and the corresponding finite difference 
(fd) approximations used in practical calculations that we want to understand 
the mechanism responsible for the suppression of van der Waals loops and the 
emergence of a singular limit of ht in an extended part of the phase diagram. 
A second reason for our investigation lies in our earlier work on hrt and its 
numerical side [12, 13, 14]: For all the merits of the theory, its practical appli- 
cation has been found to be troubled with two major difficulties that have been 
traced to the customary way of incorporating the core condition of vanishing 
pair distribution function for hard core reference systems, and to the numeri- 
cal properties of the equations for high compressibility, respectively. The latter 
clearly is an issue of prime importance when focussing on phase separation and 
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the immediate vieinity of the eritical point where kt diverges, and its severity 
can only be assessed on the basis of a thorough understanding of the numerical 
process in relation to the properties of the HRT PDE. Evidently, such an under- 
standing is also highly relevant to the interpretation of numerical results and 
the extraction of meaningful and reliable information from them. 

In order to shed some light on these questions, in the present report we 
study the analytical properties of the pde in the presence of a singular limit 
of Kt- Relegating some details to the appendix, after a short introduction to 
HRT itself and its conceptual basis in section 2 we present the pde, identify 
a quantity convenient for following the build-up of infinite compressibility, and 
infer the asymptotic scaling relations we base our work on. Employing a suitably 
formalized notion of smoothness, in sections 3 through 6 we find a total of 
three scenarios for the gradual build-up of infinite Ky, clarify their respective 
preconditions, and infer some of their properties with a view to an eventual 
implementation by FD methods. According to their most prominent traits we 
refer to the classes of asymptotic solutions so found as "(genuinely) smooth" 
(section 3), "monotonous" (section 5), and "stiff" or only "effectively smooth" 
(section 6), respectively. Of these, only the first seems to be compatible with 
an earlier analysis of hrt's scaling limit at first sight [6]; however, a closer 
investigation into the assumptions implicit in the simplifications made there 
leads to a reappraisal of those results that therefore cannot invalidate either of 
the remaining two candidate types of solution (section 7) . 

As the considerations outlined above do not take into account the initial and 
boundary conditions imposed on the PDE but rather concern themselves with a 
summary analysis of the asymptotic behavior of the various terms in the pde 
and of the range of solution types compatible with these, the all-important ques- 
tion of which of the scenarios captures the true behavior cannot be answered in 
the present contribution. In part II of our investigation [15] we present a host of 
numerical evidence strongly suggestive of a pde asymptotically turning stiff for 
thermodynamic states of infinite compressibility while FD methods on practi- 
cal discretization grids bring about a regularization that leads to an artificially 
smoothed solution. Assertion of stiffness also paves the way for a qualitative 
understanding of the relation between specific features of the Fourier transform 
of the potential and its numerical properties, clarifies the special standing of 
the hard-core Yukawa potential, and leads to a detailed and self-consistent per- 
ception of the process of solving the equations throughout the domain of the 
PDE. Needless to say, these findings are highly relevant for the interpretation 
of numerical results and the methods of data analysis to be applied to them, 
especially when working with non- uniform high-resolution discretization grids. 

2 Basic relations 

As a starting point, let us shortly review the concepts underlying HRT when 

applied to simple one-component fluids, recalling some of its central notions 
and establishing the equations we will base our work on; for a more detailed 
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account of the derivation, its physical justification and relation to both the IE 
formalism and RG theory as well as the modifications necessary for dealing with 
other physical systems, most notably spin models and fluid mixtures, we re- 
fer the reader to rcfs. [1, 12] and further references therein. For consistency 
with our earlier work on hrt we employ a number of notational conventions 
summarized in the appendix: Most importantly, a tilde indicates Fourier trans- 
formation, and once a symbol has been introduced we generally omit obvious 
generic function arguments. The appendix also serves as a repository for some of 
the more cumbersome analytical expressions as well as for auxiliary definitions 
and relations tangential to our reasoning but necessary to make our presentation 
self-contained. 

Working in the grand-canonical ensemble we consider a system of parti- 
cles interacting via pair-wise additive forces taken to derive from a potential 
v{r) = 11™' (r) -|- w{r). For the sake of simplicity, the potentials are assumed 
Q independent, and we restrict the reference fluid corresponding to v'^^^ alone 
to a system of hard spheres of diameter cr, i. e., w''®^(r) is infinite for r < a 
and vanishes otherwise. The perturbation w{r) and temperature T enter the 
calculation only in the combination (l){r) = —Piv{r), where /3 = l/ksT and ks 
is Boltzmann's constant. 

Based on this splitting oi v, & momentum space cutoff Q is introduced by 
the device of a rather artificial [16] Q dependent potential v'^^\r) obtained from 
v{r) by the elimination from w of all Fourier components w{k) with k < Q, i. e., 



Clearly, the reference and target systems with potentials v^'^ and v are obtained 
in the limits Q oo and Q ^ 0, respectively. A rather intricate analysis of a 
resummed perturbation expansion for the properties of the (Q — AQ) system 
in terms of those of the Q system at the same temperature and density in the 
limit AQ — » finally yields a non-terminating hierarchy of first-order ordinary 
differential equations (odes) in Q for the free energy A'-'^^g) and the n particle 
direct correlation functions [3]. Formally, these differential equations allow one 
to follow the evolution of structure and thermodynamics of the Q systems when 
fluctuations of ever increasing wavelength 1/Q are taken into account, i. e., 
when Q goes from infinity to zero and transformed from t;""®^ into v. 

Such an infinite set of coupled odes is, of course, hardly tractable numer- 
ically, let alone analytically. As a remedy, a closure on the two-particle level 
resembling LOGA/orpa, eq. (A4) in the appendix, is customarily adopted and 
combined with only the first HRT equation giving the Q dependence of A^^^ . As 
demonstrated in rcf. [12], if hrt's ability to describe phase coexistence is not to 
be lost, it is vital to also incorporate a condition of thermodynamic consistency 
into the closure: In HRT's standard formulation this takes the form of the com- 
pressibility sum rule (A5) relating k^'' as obtained by differentiation of the free 
energy to the volume integral of the direct correlation function at arbitrary Q. 



V' 



= V' 



,rcf 



f(r) +wW)(r), 
(fc) : k>Q 
: k<Q. 
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Due to the density derivatives so introduced the odes at fixed g give way to a 
single PDE in Q and g for the free energy that is to be solved on the semi-infinite 
strip T> where g^in < Q < Qma.^ A oo > Q > 0. The precise choice of the initial 
and boundary conditions that remain to be imposed aX Q = oo, s.t g = g-a-im, 
and &i g = ^>max is of no importance for the remainder of this work, and we 
refer the reader to refs. [12, 13, 14] for a more detailed discussion of this point. 

Discretization of the pde for the free energy so obtained is straightforward 
and yields a computational scheme that can, indeed, successfully be used for 
T > Tc, for close-to-critical and subcritical temperatures, however, attempts 
at a direct solution invariably fail to produce any results [5, 14]. In order to 
remedy this situation, Tau et al. [7] proposed an alternative formulation in 
terms of an auxiliary quantity f{Q, g) that is essentially the first Q derivative 
of the free energy; just as in our previous work on HRT [12, 13, 14], in the 
present contribution we rely on a slightly different definition for / detailed in the 
appendix, cf. eq. (A3). Further specializing to density- independent potentials 
and not explicitly including the core condition that is not expected to be relevant 
to the subject of our study, the PDE can be written in quasilinear form, 

df d^f 

g^=doo[f;Q,Q]+do2[f;Q,Q]-g^, (1) 

with initial and boundary conditions that directly follow from those imposed in 
the original formulation. 

The rather lengthy expressions for the coefficients rfoo and (io2 of eq. (1) are to 
be found in the appendix, as are the defining relations for a number of auxiliaries. 
Among these the quantity e(Q, g) = s{Q, £>) -|- 1 is of particular relevance to our 
reasoning: Essentially the exponential of /, it turns out proportional to the 
isothermal compressibility of the fully interacting system, cf. eq. (A6). Infinite 
kt therefore directly implies attendant divergences at Q = in e, e, and /, and 
we have to study the large-e behavior of the pde if we are to understand the 
description of the critical region afforded by hrt on the level of the pde. 

On the other hand, f{Q, g) is guaranteed by the construction of the hrt 
hierarchy to be continuous and finite for every non-vanishing cutoff Q and to 
coincide with its limit from above at Q ^ wherever that limit exists. In other 
words, for thermodynamic states (T, g) within the coexistence part of the phase 
diagram / must take on large finite values for sufficiently small non- vanishing Q, 
and it must diverge for Q ^ 0: If kt{8) is infinite, for every threshold F there 
is a corresponding cutoff Qf{q) > such that f{Q,g) > F for all Q < Qf{q), 
or, less formally, 

{limg^o /(Q, g) =CG 
f{Q, g) < oo : Qa > (2) 
/(0,£>)>l:Qc7<l, 

which is the very basis of the main arguments of the present report. Considering 
large e{Q, g) and assuming / to diverge more strongly than just logarithmically 
{v. L, section 5), inspection of the expressions in the appendix shows both of 
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the cocfScicnts of the quasihnear PDE (1) to be of first order in e whereas / is 
essentially the logarithm of e: 



do2 = 0(e), 
rfoo = 0(e), 

/ =o(i). 



(3) 



O and its companion o are the usual Landau symbols, and here as in the 

remainder of the present series of reports x = Oiy"") is taken to actually mean 
that a is the infimum of all h for which x = o(y^) when the implied limit is j/ ^ 
00, or else the supremum when considering y — > 0. Furthermore, qualification 
of some quantity x as "essentially" independent of some other quantity y is 
equivalent to the characterization of x as of order 0(1) in y which does not 
rule out a weak, say, logarithmic y dependence of x. Unless explicitly stated 
otherwise, all orders cited are in terms of e. 

In order to understand the properties of the solution of the pde around 
some point {Q, g) G T> we still need to supplement eq. (3) with an analogous 
characterization of the behavior of the remaining term on the right hand side 
of eq. (1), viz., d^f/dg^. Lacking any a priori information to guide us, in this 
report we adopt the simple ansatz 



a choice that is sufficiently general to admit a consistent description both of 

the behavior of the exact solution and of the computational process yielding 
a numerical approximation of it. Inserting eqs. (3) and (4) into the pde (1), 
comparison of both sides of the equation immediately yields 



only for r = is there the added possibility of a cancellation of the leading 
terms on the right hand side of eq. (1) that might give rise to an s less than 
unity, including s — r — 0. At any rate, neither r nor s may be negative. 

Equipped with eqs. (3) through (5) we are now in a position to gain a 
better understanding of the pde's workings. Although the analytical expressions 
presented and orders cited below hinge on the assumption (4), most of our 
conclusions are expected to remain qualitatively valid even in more general 
situations, c/. section V in part II [15]. 

3 Genuinely smooth solution 

From the pde itself the Q and g scales characteristic of the variations of /(Q, g) 
in the integration domain V are not apparent: In particular, we cannot say a 
priori whether they remain bounded from below in an essentially e and, hence. 




r > 



(4) 




(5) 



s ~ max(l, 1 + r) = 1 + r; 
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Figure 1: Sketch of the auxiliary function /(Q, g) for fixed and sufficiently small 
Q in the genuinely smooth scenario: ^imin and Qmax are the densities where 
boundary conditions must be imposed upon the solution of the PDE; the density 
range where / is large extends from to £>2; and there is a single, rather flat 
maximum oi f ioi Qi < g < Q2- 

temperature independent way throughout V or else scale like some inverse power 
at least of e and so become arbitrarily small in part of V whenever the temper- 
ature falls below Tc- Restricting ourselves to analytical considerations, in this 
report we will study and present arguments for both of these types of solution, 
referring to them as smooth and non-smooth, respectively. Clearly, this dis- 
tinction is highly relevant to the numerics and therefore of immediate practical 
interest: After all, smoothness implies that finite difference (fd) schemes are in 
principle well applicable to the pde at hand whereas otherwise local truncation 
errors will generally be Tinboundcd and at least an estimate of the global error 
incurred must be obtained a posteriori in order to gauge the significance of any 
information extracted from FD approximations of the PDE. 

For this section turning to the assumption of / being smooth in the sense 
stated — an assumption that we will repeatedly refer to as leading to the "(gen- 
uinely) smooth" scenario — , we immediately conclude that r = s = 0: If the 
Q and g scales set by / are bounded from below in the manner indicated, for 
any linear differential operator L wc can choose e independent step sizes AQ 
and Ag such that an fd approximation of Lf becomes accurate throughout P. 
As the estimate so obtained is just a linear combination of / values sampled 
few AQ and Ag apart, Lf is bound to scale like / and is thus of order 0(1). 
Specializing to L = d^/dg^ and to L = d/dQ we immediately conclude that 
r = and s = 0, respectively. — As a corollary we note that growth of / in 
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proportion to an inverse power of the cutoff, / oc at fixed q witli a > 0, 

always falls into this class of solutions. 

The main advantage afforded by the presupposition of smoothness is that 
it allows us to understand on the level of the pde the build-up of infinite / 
and kt close to the critical point and in the coexistence region of the phase 
diagram: Let us assume that, at some fixed and sufficiently small cutoff Q, 
f{Q,g) is a function of g like that sketched in fig. 1: continuously differentiable 
for all densities, convex and large for densities in a range with approximate 
boundaries q\ and §2 but rather small elsewhere. In the region of large / and e 
where asymptotic reasoning along the lines of eqs. (3) through (5) is applicable, 
the PDE coefficient doo is dominated by terms related to the Q dependence of the 
Fourier transforms of the potential and of the reference system direct correlation 
function, 

cf. eq. (Al). doo is therefore expected to be negative for most cutoffs less than 
the position of the first minimum of 0, and certainly for very small Q: It is by 
this standard that the cutoff under consideration must be "sufficiently small" 
as stated before. — Trivially, we also see from eq. (Al) that rfo2 is non-positive 
throughout 2?, and negative for e 0, as i^o must be positive for the PDE to be 
stable at all, cf. section 2.4.1 of ref. [14] as well as refs. [1, 7]. 

As for the remaining term on the right hand side of eq. (1), d^f/dg^ is 
negative throughout most of the density interval of large /, and positive on 
either side as / falls to small values. Furthermore, the maximum of / displayed 
in fig. 1 is rather flat so that \d'^f/dg'^\ is quite small there. Assuming this 
curvature so small that the do2 term in eq. (1) does not reverse the sign of the 
right hand side, we immediately conclude that df /dQ < where / is large, 
corresponding to further growth of / as Q = is approached from above. In 
this way we can understand how the PDE implements the gradual build-up of 
infinite kt and the attendant suppression of van der Waals loops in the free 
energy. On the other hand, it becomes clear that this mechanism depends on 
the preservation of the general outline of /: In particular, the stable growth 
of / is likely to break down once flatness is lost in the central part of the 
interval [gi,g2]- A closer look at the do2 term immediately reveals that it acts 
to stabilize this feature of the form of /: Taken by itself, rfoo strongly favors / 
to grow most rapidly close to its maximum. This, however, is prevented by the 
considerable increase in the negative density curvature of the solution it would 
entail: The correspondingly large contribution do2 {d'^f/dg'^) > to df /dQ 
effectively counteracts doo and so ensures that f{Q,g), gi < g ^ g2, remains 
flat even as it grows, just as postulated at the outset. Close to g^ and g2, on 
the other hand, the role of the rfo2 term is quite different: There the curvature 
d^f /dg^ turns positive so that both terms on the right hand side of eq. (1) 
now contribute to the growth of /, thereby rendering the transition to small / 
ever more sharply defined as Q progresses towards zero. By the same token, 
intermediate minima in the density range of large / are also dissolved by the 
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do2 term, which relaxes the precondition of conformance with fig. 1 somewhat. 

This intuitively appearling concept of a stable mechanism of growth of /, as 
we will henceforth refer to it, so allows us to understand on the level of the pde 
the emergence of infinite kj^ within a clearly demarcated density interval [qv^ 
In addition, as the coexisting densities Qy and qi are obtained as the limits of 
Qi and ^2 for Q ^ 0, anything but continuous inverse compressibility is hard to 
accomodate within this picture, which agrees well with the known coincidence 
of the HRT binodal and spinodal for space dimensionality d = 3 [6]. The key 
role played by the do2 term in stabilizing the solution and locating the densities 
of the pure phases once more highlights the importance of the thermodynamic 
consistency condition (A5) underlying the transition from odes at fixed g to a 
PDE over all of V: Indeed, it should come as no surprise that application of an 
approximation incompatible with the afore-mentioned compressibility sum rule 
yields pathological results including negative compressibility [12]. 

Without evidence to the contrary it is then tempting to subscribe to the 
seemingly natural assumption of a smooth solution, r = s = 0: Not only does 
this provide the basis for understanding the most salient features of hrt as just 
discussed, it is also what has been found in an early analysis of hrt's scaling 
limit [6], cf. section 7 below. Unfortunately, however, this view is not entirely 
unproblematic as r and s must both vanish, v. s.: As stated in connection with 
eq. (5), this directly implies that the right hand side of eq. (1) is affected by 
massive cancellation so that the sum of two terms of order 0(e) each is reduced 
to order 0(1). Indeed, the relevant signs — do2 < throughout V [eq. (Al)], 
doo < for small Q [eq. (6)], and d^f/dg'^ < for most densities in [£>i,p2] 
(fig. 1) — do allow such a cancellation. On the other hand, the rapid growth of 
the coefficient functions doi = 0(e) places rather stringent constraints on the 
shape of / at constant cutoff: According to eqs. (1) and (3) to (5), / = 0(1) 
must deviate by terms of order 0{l/s) at most from an exact solution f{Q,g) 
of the non-linear ode 

^ doo[f;Q,g] 
V do2[f;Q,g] 

for gi < g < g2. Let us now suppose that the initial and boimdary conditions 
imposed on the hrt pde actually lead to such a near-solution of the above ODE 
at some sufficiently small Q, assessment of the likelihood of which falls outside 
the scope of this report. In this case consistency with s = requires the residue 
/ — / of eq. (7) to vanish as 1/e when Q further progresses towards zero, but 
we have not been able to support this behavior of the solution on the basis 
of the explicit expressions (Al) for the doi and their properties. Whereas the 
cancellation necessary for the genuinely smooth scenario thus certainly cannot 
be ruled out, it poses far stricter preconditions than the stable mechanism of 
growth discussed before, and both its genesis and stability remain unclear at 
this point. 
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4 Refined analysis 



Looking back at the preceding section, we see that the strong point of genuine 
smoothness, viz., the stable mechanism of growth with its consequences, does not 
expHcitly depend on s = 0, the very condition that is the source of the conceptual 
difficulties just discussed. It thus seems pertinent to attempt to salvage these 
features of the smooth scenario in a less restricted setting by eliminating the 
assumption of vanishing s. But rather than separately discussing the case of 
r = 0, s = 1 we now turn to the full generality afforded by eqs. (4) and (5). 
Of course, for d^f/dg^ to be of order 0(£^) with non-vanishing r, / must be a 
rapidly oscillating function of g and we can only hope to adapt our findings if 
/ resembles fig. 1 when averaged over oscillations. 

Let us consider once more the evolution of /(Q, g) within the region of large 
£ as Q progresses towards zero: Combining eq. (5) with the definition (A2) of e 
we find 

^ = O(e-) = 0(e/^'), 



dQ 



su, 



,-.2 



0- 



Here the modified exponent s' takes into account the deviation of uo{Q) oc (fi{Q) 
from unity; for Q ^ 0, s and s' coincide by virtue of the normalization condition 
Wo(0) = 1. We now define an auxiliary quantity 

, , c?oo[/; Q, g] + do2[f; 0, g] (a^/V) 
doiQ, Q) = (8) 

= 0(1). 

Existence of do is guaranteed, and the pde (1) is fully equivalent to 

|§ = 4e/% (9) 

which no longer involves a derivative with respect to g. The signs of do and 
df /dQ always coincide. 

If only for the moment we assume s' constant, eq. (9) can be formally inte- 
grated to yield the solution at all Q given the initial condition that / be /i at 
cutoff Q\, viz., 



f{Q, q) = In (^e-f^^' + s' f^'doiq, g) dgj , 



(10) 



which is valid only if s' a s^^ does not vanish, i. e., if s > and 4> ^- The 

latter condition effectively restricts our arguments to cutoffs somewhat below 
^, the smallest of the zeros Q^^ {i = 1,2, . . ., Q^- < ^_|_^) of the Fourier 

transform ^ of the potential; in the limit Q — > that we are interested in this 
certainly poses no problem. With this proviso eq. (2) can be re-cast in form of 
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the relations 



/ do{q) dq < — — for < Qi < Q, 

Jfs' (11) 

do(q) dq = — ; — for kt = oo; 

s' 

for the thermodynamic states of interest both sides of the above inequahty are 
positive if Q — Qi is comparable to Q. 

The remainder of our analytical considerations on the nature of the PDE 
for low cutoff in the critical region will be based on eq. (11). As e"^" / s' is 
strictly monotonous in s', interval arithmetic is trivial and allows us to tackle 
in a straightforward way the problem of the Q dependence of s': Given the e 
independent and rather slow variation of s'/s oc as a function of Q, for any 
cutoff interval considered the range of s' values is easily found and translated 
into an interval of e~^^ / s' . Arguments based upon eq. (11) like those we will 
present shortly are thus easily modified to take into account the non-constancy 
of s' simply by applying the least restrictive bound for any of the s' values in 
the Q interval under consideration, an operation taken for granted throughout 
the remainder of this series of reports. 




5 Monotonous growth and logarithmic singular- 
ity 

According to its definition (8), the integrand do on the left hand side of eq. (11) 
is always of order 0(1). On the other hand, for / to be a monotonous function of 
Q in the asymptotic region, its reduced slope do must always be negative there 
so that — jQ^do{q) dq = 0(1) {Q — Qi). As Qi is an essentially free parameter 
to be chosen from (0, Q], eq. (11) shows a quantity scaling like Q to be bounded 
from above by another one of order 0(e~^). i. e., 0{Q) < 0{s~^). This prompts 
us to consider a power law relation between Q and e, say, 

Q = o(£-*), 

t > s, 

corresponding to only a logarithmic singularity of / at Q = 0. In this case, 
eq. (3) no longer holds due to the cutoff dependence of the prefactors in eq. (Al): 
To leading order we find 

do2 = 0(e-i-2*), 

doo = 0(£l-*), 

/ =o(i) 

instead, and the balance equation is modified to 

0(£«) = 0(ei-*) + 0(£i-2*)0(£^). (12) 
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If the leading terms on the right hand side of eq. (12) do not cancel, eq. (5) 
is thus to be replaced by the relation s = max(l — t,l — 2t + r) with r > 
and t > s > 0. For r > t, however, the assumption of monotonous growth 
of / leads to inconsistencies: In this case, the right hand side of the pde (1) 
is asymptotically dominated by the do2 term. As rapid oscillations in the g 
direction at fixed Q provide the only way for d'^f/dg'^ to become exponentially 
large compared to / itself, the second density derivative of / is bound to be 
oscillatory in the density range gi < g < g2 just as well. Unlike /, however, 
9^f/9g^ must change its sign at every swing, which immediately carries over to 
df /dQ and do due to the asymptotic dominance of the do2 term, contrary to 
the assumption of do{q) being negative for all q < Q. 

We are thus left with the possibility of r < t, s = 1 —t; combining this with 
the conditions r > 0, s > 0, and t > s and eliminating t we find the admissible 
exponent ranges 

<r< 1, 

<s< min(i,l -r) < i. 

For r = t, there is the added possibility of cancellation of the leading terms on 
the right hand side of eq. (12); again eliminating t from the relations < r = t, 
< s < 1 — t, and t>swe obtain 

<r< 1, 

<s< min(r, 1-r) < i. 

At any rate, non- vanishing exponent s is compatible with monotonous growth 

of / only in the case of a merely logarithmic singularity at Q = 0, and in this 
case there is an upper bound of ^ for s that holds irrespective of whether the 
leading terms on the right hand side of eq. (12) cancel. Furthermore, in this 
monotonous growth scenario 0{Q) < 0(e~^), v. s., so that finally s must 
tend to a finite, possibly vanishing limit for Q ^ 0. 

6 Effective smoothness from stiffness 

An alternative is to give up the monotonicity assumption and allow do to alter- 
nate in sign: This opens up the possibility of a partial cancellation of the positive 
and negative contributions to the integral in eq. (11), and the average of do no 
longer has to be of order 0(1) even though do itself still is. As we will now show, 
this situation implies that the HRT PDE turns stiff for thermodynamic states of 
infinite compressibility, with r > and s > 1. In numerical applications, on 
the other hand, discretization of the pde on practical grids by necessity induces 
an artificial smoothing of the numerical solution that is then characterized by 
vanishing effective exponents refj and Seff provided the computation is able to 
reach the limit Q ^ at all. 

Stiffness of the pde in that part of T> where the divergence of kt is built up 
follows from cq. (11): In the scenario outlined above, the bound e~-^^ /Qs' = 
0{e~^)/Q on the mean of —do over the interval (0,(3] implies that do is a 
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rapidly oscillating function of Q both amplitude and period of which scale like 
l/e^Q = 0(£~*) [the possibility of a non- monotonous logarithmic divergence, 
Q ~ e~* with < t < s, may be accounted for by replacing s with s — t > 
0]. Considering cq. (9), for e ^ / ^ 1 the solution is then characterized by 
oscillations with similar scaling properties superimposed upon a large, smooth, 
and most likely monotonously growing regular part. These oscillations in the Q 
direction are naturally accompanied by oscillations in the g direction on density 
scales of order 0{e~^^'^), r > 0: Vanishing r would require the evolution of / in 
Q at densities at fixed separation to remain synchronous over many, viz., on the 
order of Q undulations of the function; this seems highly unlikely, and r = 
must be unstable in the setting under discussion. From eq. (5) we then obtain 
r > and s = r + 1 > 1, where the possibility of a weak Q and g dependence 
of r and s must also be anticipated. 

The scenario just discussed differs from genuine smoothness and monotonous 
growth in two important ways: Firstly, the singularity of / may be stronger than 
before as there is now no need for the cancellation required in section 3, nor 
does boundedness of the mean of do imply boundedness of the mean slope of 
/ as in section 5. Secondly, non- vanishing exponents r and s mean that the 
scales characteristic of the variation of / and, hence, the grid spacings AQ 
and Ag appropriate for the discretization of the pde become arbitrarily small 
as diverging isothermal compressibility is built up, which obviously has grave 
repercussions for the applicability of fd methods to the pde at hand. 

This last point is worth considering in some detail: Clearly, for the fd equa- 
tions to be a good approximation of the pde the step sizes must go to zero as 
appropriate inverse powers of e determined by r, s, and the orders of the local 
truncation error in AQ and Ag. Excluding a very weak divergence of /, however, 
the rapid growth of e certainly renders so fine a grid impractical [12], and below 
some cutoff Qaq {T, g) > the numerics can no longer follow the evolution in 
Q of the true solution of the pde, nor that in g below some QAe{T,g) > 0. 
Except right at gi and £>2, both Qax siie rather well defined because of the 
non-zero exponents s > r > 0. Below the QAx^ the basic assumption underly- 
ing any fd method, viz., applicability of Taylor expansions of rather low order, 
tantamount to smoothness of the solution on the scales set by the discretization 
grid spacings, no longer holds so that a discretization derived on this basis is 
used outside its range of validity. A chosen implementation of the theory may 
therefore fail to proceed to Q < Qax altogether due to inappropriately large 
step sizes: In this case no solution is ever obtained at Q = below the critical 
temperature and neither criticality nor phase separation can be described, just 
as reported in ref. [5] and appendix B.l of ref. [14]. Alternatively, judicious 
choice of the formulation of the theory and of its discretization may allow the 
numerical scheme to solve the fd equations all the way to vanishing cutoff; it 
is this very property of eq. (1) that prompted adoption of the quasilinear form 
of the HRT PDE in the first place, cf. section 2. Taylor expansion arguments no 
longer being applicable, however, a priori bounds on the local truncation error 
are not available, and the global error may well be substantial. Furthermore, as 
any solution generated in a FD calculation on a given grid is well represented on 
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that same grid by definition, wc can use the same arguments as in section 3 to 
show that the numerical results obtained with step sizes of order 0(1) necessar- 
ily reproduce vanishing effective exponents r^s = s^s = 0. This reduction of the 
exponents from s > r > for the exact solution corresponds to a smoothing of 
/ that presumably weakens the singularity, v. s., and also suppresses any oscil- 
lations on cutoff and density scales smaller than the step sizes. Such smoothing 
is very well known in the numerics of PDEs and in fact forms the conceptual 
basis for the highly efHcient multi-grid methods for the solution of integral and 
differential equations [17, 18]. 

In this "stiff," or only "effectively smooth" scenario, any results obtained 
by FD methods then necessarily realize the genuinely smooth case of section 3, 
admitting, i. a., a numerical solution / oc 1/Q. In this case the artificial smooth- 
ing attendant to the reduction of the exponents directly affects the solution for 
T <Tc at any cutoff lower than (5smooth(^r, g) = max((5AQ, Qaq), and the final 
results for Qy < Q < Qi- Given the do2 term of the pde, however, such drastic 
qualitative (rofr ^ r) and quantitative {scs ^ s) changes inside the region of 
large e must have an effect on the location of the binodal as well as on the solu- 
tion for the pure phases, too. Restriction of large / to only a rather well-defined 
density interval [pi,P2]j the pde's parabolic character, and the boundary con- 
ditions imposed at Q^in and Pmax nevertheless inspire some hope that the error 
incurred outside the binodal might be limited. Indeed, in the presence of stiff- 
ness and smoothing this very hope must underlie any attempt at a numerical 
solution of the hrt equations and certainly has to be justified at least a poste- 
riori by combining related calculations and testing the internal consistency of 
the results obtained [12, 13, 14]. 



7 Stiffness and the scaling limit of HRT 

So far our considerations have led us to the formulation of three different pos- 
sibilities for the asymptotic behavior of the hrt pde for high compressibility 
states characterized by true smoothness of the solution, by monotonous growth 
leading to only a logarithmic singularity, and by stiffness giving rise to effective 
smoothness of /, respectively. 

At first sight, however, anything but the genuinely smooth solution, r = s = 
0, seems to be at variance with the detailed analysis of hrt's scaling limit in 
ref. [6] : In section III B of that report we find the prediction that the quantity 
denoted uq there (corresponding to — /) should be quadratic in the density-like 
variable x and scale like Q^^"^^^^ in d dimensions; in our notation, / should thus 
grow like 1/Q in d = 3 dimensions, corresponding to r = s = 0, cf. section 3. 

In order to resolve this seeming contradiction let us have a closer look at the 
reasoning of ref. [6] for T < T^. The pivotal relation is eq. (3.9) there, viz., 

For this equation, the direct analogue of our eq. (1), the authors of ref. [6] invoke 
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the validity of "neglecting exponentially small terms" to justify replacing its 
left hand side by nought; the asymptotic solution proportional to Q"*^''"^' then 
follows immediately by separation of the cutoff and density dependencies. But 
this replacement is legal only if the slope of uq actually remains small in modulus 
compared to its exponential or, translating to our notation, if (1/e) {df /dQ) 
tends to zero as / goes to +00, i. e., if s < 1 as is the case for the solution cited 
where it actually vanishes. Conversely, cq. (13) can also be read to indicate 
that its left hand side will be small if and only if the pde's initial and boundary 
conditions have caused uq at the cutoff considered to almost exactly solve a 
simple ODE in the density-like variable x, viz., 

which is completely equivalent to the cancellation requirement of section 3, 
eq. (7) in particular. The above directly shows that d'^ug/dx'^ is proportional 
to Q^~'^ oc Uq and so of order 0(1) in e~"'?; in terms of / this means that 
d'^f/dg^ is of order 0(1) in e and, hence, that r = 0. 

As the assumptions both of smoothness in the g direction, r = 0, and of 
cancellation of the leading terms on the right hand side of eq. (1), s < 1, are 
thus built right into the derivation of the scaling solution of ref . [6] , conformance 
of the latter with the smooth scenario is hardly surprising and certainly does not 
rule out any of the other possibilities considered in the present report. As for the 
numerics, both genuine (r = s = 0) and effective (reff = Sgg = 0) smoothness 
suggest an / that asymptotically grows like 1/Q, although in the latter case 
this docs not reflect the properties of the exact solution of the pde. We are 
thus led to a reappraisal of the analysis presented in section IIIB of ref. [6] 
as certainly applicable to the fdes solved numerically but not to the true pde 
unless r = s = 0, i. e., unless the cancellation requirement of section 3 is 
also met; with this proviso the conclusions of ref. [6] remain largely unaffected. 
As for the intuitively appealing stable mechanism of growth responsible for 
suppression of van der Waals loops, it now applies to all three scenarios, if only 
for the numerical process unless r = 0; further support for it derives from the 
^-dependence of uq that conforms with the sketch of fig. 1 whereas eq. (7) might 
well be compatible with a more general outline. 

In view of this resolution of the seeming contradiction between the earlier 
analysis of hrt's scaling limit and the possibility of non-zero exponents r and s 
none of the three scenarios can be excluded from further consideration at this 
point. On the other hand, the rather summary analytical considerations pre- 
sented here can only identify the solution types consistent with the asymptotic 
properties of the pde for infinite compressibility and explore their preconditions 
and consequences, which has been the subject matter of the present report. The 
all-important question of which of them is in fact realized in actual applications 
of the theory, however, crucially depends on global aspects of eq. (1) in all of V 
and certainly cannot be answered without considering the influence of the initial 
and boundary conditions that, after all, uniquely determine f{Q, g) throughout 
the integration domain. To arrive at a decision we therefore need to solve a 
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discretized version of the pde, scrutinize the computational process, and inter- 
pret the numerical evidence so gained, a task we will undertake in part II of our 
investigation [15]. 
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A Notational conventions and additional rela- 
tions 

A most detailed account of the process of re- writing the pde for the free energy 
into the quasilinear eq. (1) can be found, alongside the explicit expressions for 
the PDE coefficients doi, in appendix A of ref. [14]. Further specializing these 
to take into account the elimination of the core condition and the purported 
density-independence of the potential we obtain 



doo = +7^ 




2/ 



do2 



(Al) 



In the above expressions we have suppressed the obvious function arguments Q 
and Q, and we rely on the notational convention introduced in our earlier work 
on HRT according to which superscripts indicate the system a quantity refers to 
and a tilde signals Fourier transformation. We also make use of the following 
auxiliary quantities: 

wo(r-) = (t){r)/4>o, 
4>o = ^(0), 

}C{k,Q) = -}-+cf{k,g), (A2) 

IneiQ, g) = f{Q, g) ul{Q) - 4>{Q)/t{Q, q), 
£-(g, g) = £(Q, g) - 1. 

As far as the definition of K, is concerned, the ideal gas term involving —1/g 
is customarily included in the definition of the hard sphere reference system 
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direct correlation function d^^ throughout much of the htcrature on HRT. Fur- 
thermore, when exphcitly taking into account the core condition /C is typically 
augmented by a truncated series expansion of the deviation of the direct cor- 
relation function of the Q system inside the core; from that of the reference 
system, a scheme designed to allow one to side-step the need for costly numeri- 
cal Fourier transformations that would otherwise arise from the inversion of the 
Ornstcin-Zcrnikc equation [5, 12]. 

The connection to the thermodynamics as well as to the original formulation 
of the theory is afforded by the relations linking / and the auxiliary quantities 
defined in cq. (A2) to the first derivative with respect to Q of the free energy 
A^'^^ of the system at cutoff Q: From the expressions given in ref. [12] one may 
easily show that the evolution of A^*^^ is given by 

d pA^Q^ Q 



= (ins - 4> gj for Q > 0, 



dQ V 



(A3) 

'0 



V V 2 

The structure of the Q system follows from the closure relation that imposes a 
direct correlation function C2^^ (r, g) of the form 

JQ) _ ^ref , - (Q) , (Q) 



+ ir' uo (A4) 

parameterized by the single scalar Jq^\q) that is determined throughout V 
from thermodynamic consistency in the form of the compressibility sum rule 

In the spirit of LOGA/orpa, when implementing the core condition eq. (A4) 
is taken to hold only for r > a, and the direct correlation function inside the 
core is optimized so as to approximately minimize the pair distribution function 
there. As shown in appendix A. 3 of ref. [14], the isothermal compressibility 
of the Q system is also readily evaluated and, in the limit Q 0, found to 
reduce to 
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